### Relevant packages----------
library(haven)
library(car)
library(dplyr)
library(glmnet)
library(lmtest)
library(sandwich)
library(foreign)
library(multcomp)
library(sciplot)
library(ggplot2)
library(coefplot)
library(texreg) 

### Import relevant data--------------

##set your working directory to the location of these files

## Pretest dataset
caliphatedatfinal <- read_stata("caliphate-knowledge.dta") 

## Journalistic corrections dataset
datfinal <- read_stata("journalistic-corrections-final.dta")

### Baseline knowledge test (pre-test)------------
## Testing for recognition of the Caliphate podcast story

## Remove initial rows from Qualtrics data (after submission cutoff)
caliphatedatefinal <- caliphatedatfinal[-1, ][-2, ][-3, ]

## Remove 2nd row (still appearing after initial screenout)
caliphatedatafinal <- caliphatedatefinal[caliphatedatefinal$Status != 1,]

## Remove inattentive respondents
caliphatedatafinal_att <- caliphatedatafinal[(caliphatedatafinal$Q829_2 == 6 | caliphatedatafinal$Q829_2 == 7 | caliphatedatafinal$Q829_2 == 8) & (caliphatedatafinal$Q34_6 == 1 | caliphatedatafinal$Q34_6 == 2 | caliphatedatafinal$Q34_6 == 3), ]
# Drops 19 / 271 respondents (final N = 252)

## Calculate knowledge levels
## Identifying NYT issued major retraction (Q1005_4)
table(caliphatedatafinal_att$Q1005_4)
# 54 respondents attribute to NYT (were allowed to select multiple outlets)

## Identifying retraction topic (Response #1 to Q1146)
table(caliphatedatafinal_att$Q1146)
# 6 respondents identify "ISIS whistleblower" as topic of retraction

## Description of retraction
# Print only responses from those who identified NYT, ISIS whistleblower, and provided more detail (5 respondents)
print(caliphatedatafinal_att$Q1008[caliphatedatafinal_att$Q1008 != ""])
# Only 2 / 252 respondents (0.8 %) recall that the retraction had to do with an unreliable source in the NYT's Caliphate podcast

### Main study code----------------
## Screen out attention check failure
datafin <- datfinal[(datfinal$Q829_2 == 6 | datfinal$Q829_2 == 7 | datfinal$Q829_2 == 8) & (datfinal$Q830_6 == 1 | datfinal$Q830_6 == 2 | datfinal$Q830_6 == 3), ]

## Twitter feed response time test
## Response Time 1: Time participants spent viewing the initial stimulus (given to all participants)
datafin$response_time_1 <- NA
datafin$response_time_1 <- datafin$Q2101_Page_Submit

## Response Time 2: Time participants spent looking at their particular set of treatment tweets (aggregated into one table)
datafin$treatment_response_time_2 <- NA
datafin$treatment_response_time_2 <- ifelse(is.na(datafin$Q1054_Page_Submit) == FALSE, datafin$Q1054_Page_Submit, 
                                            ifelse(is.na(datafin$Q1056_Page_Submit) == FALSE, datafin$Q1056_Page_Submit, 
                                                   ifelse(is.na(datafin$Q1058_Page_Submit) == FALSE, datafin$Q1058_Page_Submit,
                                                          ifelse(is.na(datafin$Q1060_Page_Submit) == FALSE, datafin$Q1060_Page_Submit,
                                                                 ifelse(is.na(datafin$Q1062_Page_Submit) == FALSE, datafin$Q1062_Page_Submit, NA)))))
summary(datafin$response_time_1)
summary(datafin$treatment_response_time_2)

## Recoding treatment assignment
datafin$condition_num <- NA
datafin$condition_num[datafin$Condition == "NoCorrection"] <- 1
datafin$condition_num[datafin$Condition == "MediaCorrection"] <- 2
datafin$condition_num[datafin$Condition == "3Party"] <- 3
datafin$condition_num[datafin$Condition == "Media3Party"] <- 4
datafin$condition_num[datafin$Condition == "3PartyMedia"] <- 5

# Verifying random assignment
table(datafin$condition_num)

## Treatment assignment indicator variable
# Initializing as '0' rather than 'NA' to handle missing values 
# (R defaults =/= Stata defaults)
datafin$no_correction <- 0
datafin$news_correction <- 0
datafin$third_correction <- 0
datafin$news_third_correction <- 0
datafin$third_news_correction <- 0
datafin$no_correction[datafin$Condition == "NoCorrection"] <- 1
datafin$news_correction[datafin$Condition == "MediaCorrection"] <- 1
datafin$third_correction[datafin$Condition == "3Party"] <- 1
datafin$news_third_correction[datafin$Condition == "Media3Party"] <- 1
datafin$third_news_correction[datafin$Condition == "3PartyMedia"] <- 1

## Recoding pretreatment outcome variables
## Belief accuracy
# Note: scale inverted from Qualtrics survey output 
# Higher values indicate more accurate beliefs (i.e. that story was false)
names(datafin)[137] <- "belief_accuracy_pre"
datafin$belief_accuracy_pre <- ifelse(datafin$belief_accuracy_pre == 1, 4,
                                      ifelse(datafin$belief_accuracy_pre == 2, 3,
                                             ifelse(datafin$belief_accuracy_pre == 3, 2,
                                                    ifelse(datafin$belief_accuracy_pre == 4, 1, NA))))

## Trust of specific news outlet (the CBC) BEFORE initial stimulus
datafin$news_trust_prestimulus <- NA
datafin$news_trust_prestimulus[datafin$Q1127_3 == 1] <- 1
datafin$news_trust_prestimulus[datafin$Q1127_3 == 2] <- 2
datafin$news_trust_prestimulus[datafin$Q1127_3 == 3] <- 3
datafin$news_trust_prestimulus[datafin$Q1127_3 == 4] <- 4

## Trust of specific news outlet (the CBC) AFTER initial stimulus but BEFORE treatment
table(datafin$Q1139_3)
datafin$news_trust_pre <- NA
datafin$news_trust_pre[datafin$Q1139_3 == 1] <- 1
datafin$news_trust_pre[datafin$Q1139_3 == 2] <- 2
datafin$news_trust_pre[datafin$Q1139_3 == 3] <- 4
datafin$news_trust_pre[datafin$Q1139_3 == 5] <- 3

## Political Efficacy 
# Renaming individual questions
names(datafin)[120] <- "poli_efficacy_pre1"
names(datafin)[121] <- "poli_efficacy_pre2"
names(datafin)[122] <- "poli_efficacy_pre3"

# Singular metric (mean of three questions)
datafin$poli_efficacy_pre <- NA
datafin$poli_efficacy_pre <- rowMeans(datafin[120:122], na.rm = TRUE)

## Coding moderating variables
## Republican identification
datafin$Republican <- NA
datafin$Republican[datafin$Q778 == 2 | ((datafin$Q778 == 3 |datafin$Q778 == 4) & datafin$Q778 != 1)] <- 0
datafin$Republican[datafin$Q778 == 1 | datafin$Q782 == 1] <- 1

## Generalized trust in the news/mass media
datafin$media_trust <- NA
datafin$media_trust[datafin$Q1179 == 1] <- 1
datafin$media_trust[datafin$Q1179 == 2] <- 2
datafin$media_trust[datafin$Q1179 == 3] <- 3
datafin$media_trust[datafin$Q1179 == 4] <- 4

## Twitter usage
#Initialized as 0 to use as bargraph grouping factor
datafin$twitter_user <- 0
datafin$twitter_user[datafin$SocialMedia_2 == 1] <- 1

## Outcome measures
## Belief accuracy
datafin$belief_accuracy <- NA
datafin$belief_accuracy <- datafin$Q1052_3
datafin$belief_accuracy <- ifelse(datafin$belief_accuracy == 1, 4,
                                  ifelse(datafin$belief_accuracy == 2, 3,
                                         ifelse(datafin$belief_accuracy == 3, 2,
                                                ifelse(datafin$belief_accuracy == 4, 1, NA))))

## News outlet trust
datafin$news_trust <- NA
datafin$news_trust[datafin$Q1144_3 == 1] <- 1
datafin$news_trust[datafin$Q1144_3 == 2] <- 2
datafin$news_trust[datafin$Q1144_3 == 3] <- 3
datafin$news_trust[datafin$Q1144_3 == 4] <- 4

## Political efficacy
names(datafin)[215] <- "poli_efficacy_1"
names(datafin)[216] <- "poli_efficacy_2"
names(datafin)[217] <- "poli_efficacy_3"

datafin$poli_efficacy <- NA
datafin$poli_efficacy <- rowMeans(datafin[215:217], na.rm = TRUE)

## Control variables (selected for regressions by LASSO)------
## Agegroup 
# Factor variable (1/0)
datafin$age_18 <- 0
datafin$age_18[datafin$Q808 == 2] <- 1

datafin$age_25 <- 0
datafin$age_25[datafin$Q808 == 3] <- 1

datafin$age_35 <- 0
datafin$age_35[datafin$Q808 == 4] <- 1

datafin$age_45 <- 0
datafin$age_45[datafin$Q808 == 5] <- 1

datafin$age_55 <- 0
datafin$age_55[datafin$Q808 == 6] <- 1

datafin$age_65 <- 0
datafin$age_65[datafin$Q808 == 7 | datafin$Q808 == 8 | datafin$Q808 == 9] <- 1
# 65 and older grouped together

## Education (college graduate indicator - bachelor's or higher)
# Factor variable (1/0)
datafin$college_educ <- NA
datafin$college_educ[datafin$Q776 < 5] <- 0
datafin$college_educ[datafin$Q776 >= 5] <- 1

## Gender (male)
# Factor variable (1/0)
datafin$male <- 0
datafin$male[datafin$Q768 == 1] <- 1

## Political interest
# Five-point scale (higher values indicate greater interest)
# Note: Reversed from Qualtrics output
datafin$poli_interest <- NA
datafin$poli_interest[datafin$Q790 == 5] <- 1
datafin$poli_interest[datafin$Q790 == 4] <- 2
datafin$poli_interest[datafin$Q790 == 3] <- 3
datafin$poli_interest[datafin$Q790 == 2] <- 4
datafin$poli_interest[datafin$Q790 == 1] <- 5

## Race (Non-Hispanic White)
# Factor variable (1/0)
datafin$nonhispanicwhite <- NA
datafin$nonhispanicwhite[datafin$Q770_2 == 1 | datafin$Q770_3 == 1 |datafin$Q770_4 == 1 | datafin$Q770_5 == 1 |datafin$Q770_6 == 1 | datafin$Q772 == 1] <- 0
datafin$nonhispanicwhite[datafin$Q770_1 == 1 & is.na(datafin$Q770_2) == TRUE & is.na(datafin$Q770_3) == TRUE & is.na(datafin$Q770_4) == TRUE & is.na(datafin$Q770_5) == TRUE & is.na(datafin$Q770_6) == TRUE & datafin$Q772 != 1] <- 1

## Descriptive statistics------------
## Pretreatment variables
summary(datafin$belief_accuracy_pre)
summary(datafin$news_trust_prestimulus)
summary(datafin$news_trust_pre)
summary(datafin$poli_efficacy_pre)

## Outcome variables
summary(datafin$belief_accuracy)
summary(datafin$news_trust)
summary(datafin$poli_efficacy)

## Pretreatment and outcome variables (proportions)
table(datafin$belief_accuracy_pre)
table(datafin$belief_accuracy)
table(datafin$news_trust_prestimulus)
table(datafin$news_trust_pre)
table(datafin$news_trust)

## Treatment effect moderators
table(datafin$Republican)
table(datafin$media_trust)
table(datafin$twitter_user)

## Balance table stats
## Used to create Table B1
table(datafin$male, datafin$Condition) %>% prop.table(margin = 2)
mean(datafin$male)
table(datafin$college_educ, datafin$Condition) %>% prop.table(margin = 2)
mean(datafin$college_educ, na.rm=TRUE)
table(datafin$nonhispanicwhite, datafin$Condition) %>% prop.table(margin = 2)
mean(datafin$nonhispanicwhite, na.rm=TRUE)
table(datafin$age_18, datafin$Condition) %>% prop.table(margin = 2)
mean(datafin$age_18)
table(datafin$age_25, datafin$Condition) %>% prop.table(margin = 2)
mean(datafin$age_25)
table(datafin$age_35, datafin$Condition) %>% prop.table(margin = 2)
mean(datafin$age_35)
table(datafin$age_45, datafin$Condition) %>% prop.table(margin = 2)
mean(datafin$age_45)
table(datafin$age_55, datafin$Condition) %>% prop.table(margin = 2)
mean(datafin$age_55)
table(datafin$age_65, datafin$Condition) %>% prop.table(margin = 2)
mean(datafin$age_65)
table(datafin$Republican, datafin$Condition) %>% prop.table(margin = 2)
mean(datafin$Republican, na.rm=TRUE)

### ANALYSIS-------------

## LASSO for outcome measures
## Belief accuracy
# Per summary, 112 NA values in final sample (~3.9% of sample)
summary(datafin$belief_accuracy)
# Thus, using listwise deletion to only include respondents who reached outcome variable

datafinlassoFB <- datafin[is.na(datafin$belief_accuracy) == FALSE, ]
y <- datafinlassoFB$belief_accuracy
x <- data.matrix(datafinlassoFB[, c("age_18", "age_25", "age_35", "age_45", "age_55", "age_65", "college_educ", "male", "poli_interest", "nonhispanicwhite", "belief_accuracy_pre")])
# test for missing values in control variable
summary(x)
# remove missing values from control variable
x[is.na(x)] <- 0

# alpha = 1 so glmnet runs LASSO (vs. ridge, elastic net, etc.)
# using Gaussian family for continuous outcome variable scale
cv_model <- cv.glmnet(x, y, family = "gaussian", alpha=1)
best_lambda <- cv_model$lambda.min
best_modelFB <- glmnet(x, y, family = "gaussian", alpha = 1, lambda = best_lambda)
## Non-zeroed coefficients are selected by LASSO as potential controls
coef(best_modelFB)
## Relevant controls: age_65 and belief_accuracy_pre

## News trust
datafinlassoNT <- datafin[is.na(datafin$news_trust) == FALSE, ]
y <- datafinlassoNT$news_trust
x <- data.matrix(datafinlassoNT[, c("age_18", "age_25", "age_35", "age_45", "age_55", "age_65", "college_educ", "male", "poli_interest", "nonhispanicwhite", "news_trust_pre", "news_trust_prestimulus")])
summary(x) # 4 NA values
# remove NA values
x[is.na(x)] <- 0
cv_model <- cv.glmnet(x, y, family = "gaussian", alpha=1)
best_lambda <- cv_model$lambda.min
best_modelNT <- glmnet(x, y, family = "gaussian", alpha = 1, lambda = best_lambda)
coef(best_modelNT)
## Relevant controls:age_18, age_55, poli_interest, nonhispanicwhite, news_trust_pre, news_trust_prestimulus (news_trust_pre > news_trust_prestimulus)

## Political efficacy
datafinlassoPE <- datafin[is.na(datafin$news_trust) == FALSE, ]
y <- datafinlassoPE$news_trust
x <- data.matrix(datafinlassoPE[, c("age_18", "age_25", "age_35", "age_45", "age_55", "age_65", "college_educ", "male", "poli_interest", "nonhispanicwhite", "poli_efficacy_pre")])
summary(x)
x[is.na(x)] <- 0 # remove 2 NA values
cv_modelPE <- cv.glmnet(x, y, family = "gaussian", alpha=1)
best_lambda <- cv_model$lambda.min
best_modelPE <- glmnet(x, y, family = "gaussian", alpha = 1, lambda = best_lambda)
coef(best_modelPE)
## Relevant controls: age_18, age_25, age_55, age_65, college_educ, poli_interest, nonhispanicwhite, poli_efficacy_pre

## Manipulation check passage rates (overall and by condition)
# Create indicator variables
datafin$manip0 <- 0
datafin$manip0[datafin$Condition == "NoCorrection"] <- 1
datafin$manip1 <- 0
datafin$manip1[datafin$Condition == "MediaCorrection" | datafin$Group == "1"] <- 1
datafin$manip2 <- 0
datafin$manip2[datafin$Condition == "3Party" | datafin$Group == "2"] <- 1

## Strict (selected ONLY the observed treatment tweets)
datafin$manipulated <- ifelse(is.na(datafin$Q1182_1) == TRUE & is.na(datafin$Q1182_2) == TRUE & datafin$Q1182_3 == 1 & datafin$Q1182_4 == 1 & is.na(datafin$Q1182_5) == TRUE, 1, 0)
sum(datafin$manipulated,na.rm=TRUE)/length(datafin$Condition)

## No correction treatment
datafin$manipulated_t0 <- ifelse(datafin$Q1178_1 == 1 & is.na(datafin$Q1178_2) == TRUE & datafin$Q1178_3 == 1 & is.na(datafin$Q1178_4) == TRUE & is.na(datafin$Q1178_5) == TRUE, 1, 0)
sum(datafin$manipulated_t0,na.rm=TRUE)/sum(datafin$Condition=="NoCorrection",na.rm=TRUE)

## News outlet correction treatment
datafin$manipulated_t1 <- ifelse(datafin$Q1177_1 == 1 & datafin$Q1177_2 == 1 & is.na(datafin$Q1177_3) == TRUE & is.na(datafin$Q1177_4) == TRUE & is.na(datafin$Q1177_5) == TRUE, 1, 0)
sum(datafin$manipulated_t1, na.rm=TRUE)/(sum(datafin$Condition=="MediaCorrection",na.rm=TRUE)+sum(datafin$Condition=="3PartyMedia",na.rm=TRUE)+sum(datafin$Condition=="Media3Party",na.rm=TRUE))

## Third party correction treatment
datafin$manipulated_t2 <- ifelse(datafin$Q1181_1 == 1 & datafin$Q1181_2 == 1 & is.na(datafin$Q1181_3) == TRUE & is.na(datafin$Q1181_4) == TRUE & is.na(datafin$Q1181_5) == TRUE, 1, 0)
sum(datafin$manipulated_t2, na.rm=TRUE)/(sum(datafin$Condition=="3Party",na.rm=TRUE)+sum(datafin$Condition=="3PartyMedia",na.rm=TRUE)+sum(datafin$Condition=="Media3Party",na.rm=TRUE))

## HYPOTHESIS TESTING------------
## H1 (effect of treatments on belief accuracy)
# OLS regresion with robust errors using controls selected via LASSO
Treatment_BA <- lm(datafin$belief_accuracy ~ datafin$news_correction + datafin$third_correction 
                   + datafin$news_third_correction + datafin$third_news_correction
                   + datafin$age_65  + datafin$belief_accuracy_pre) 

# Adding robust errors using 'sandwich' package
coeftest(Treatment_BA, vcov = vcovHC, type = "HC1")
coefci(Treatment_BA, vcov = vcovHC, type = "HC1")

## Test for differential effects using linear combination(s) between treatments
## Both corrections (news outlet first) - news outlet correction
# Computing via glht method
# Create matrix of model coefficients to test
# (Names and positions pulled from H1 regression output)
K <- matrix(c(0,1,0,-1,0,0,0),1)
# Test generalized linear hypothesis
summary(glht(Treatment_BA, linfct = K))

# Both corrections (news outlet first) - third party correction
## Included in Table 1 (Main effects)
K <- matrix(c(0,0,1,-1,0,0,0),1)
summary(glht(Treatment_BA, linfct = K))

# Both corrections (third party first) - news outlet correction
## Included in Table 1 (Main effects)
K <- matrix(c(0,1,0,0,-1,0,0),1)
summary(glht(Treatment_BA, linfct = K))

# Both corrections (third party first) - third party correction
## Included in Table 1 (Main effects)
K <- matrix(c(0,0,1,0,-1,0,0),1)
summary(glht(Treatment_BA, linfct = K))

# Both corrections (news outlet first) - both corrections (third party first)
## Included in Table 1 (Main effects)
K <- matrix(c(0,0,0,-1,1,0,0),1)
summary(glht(Treatment_BA, linfct = K))

# Exploratory - difference in treatment effects of news_correction and third_party correction
K <- matrix(c(0,1,-1,0,0,0,0),1)
summary(glht(Treatment_BA, linfct = K))

## H2 (Effect of treatment on news outlet trust)
Treatment_NT <- lm(datafin$news_trust ~ datafin$news_correction + datafin$third_correction
                   + datafin$news_third_correction + datafin$third_news_correction
                   + datafin$age_18 + datafin$age_55
                   + datafin$poli_interest + datafin$nonhispanicwhite 
                   + datafin$news_trust_pre + datafin$news_trust_prestimulus)

coeftest(Treatment_NT, vcov = vcovHC, type = "HC1")
coefci(Treatment_NT, vcov = vcovHC, type = "HC1")

# Same set of linear combinations for news trust
# Both corrections (news outlet first) - third party correction
## Included in Table 1 (Main effects)
K <- matrix(c(0,0,-1,1,0,0,0,0,0,0,0),1)
summary(glht(Treatment_NT, linfct = K))

# Both corrections (third party first) - third party correction
## Included in Table 1 (Main effects)
K <- matrix(c(0,0,-1,0,1,0,0,0,0,0,0),1)
summary(glht(Treatment_NT, linfct = K))

# Both corrections (news outlet first) - both corrections (third party first)
## Included in Table 1 (Main effects)
K <- matrix(c(0,0,0,1,-1,0,0,0,0,0,0),1)
summary(glht(Treatment_NT, linfct = K))

# Fill out table with additional linear combinations
# Both corrections (news outlet first) - news outlet correction
## Included in Table 1 (Main effects)
K <- matrix(c(0,-1,0,1,0,0,0,0,0,0,0),1)
summary(glht(Treatment_NT, linfct = K))

# Both corrections (third party first) - news outlet correction
## Included in Table 1 (Main effects)
K <- matrix(c(0,-1,0,0,1,0,0,0,0,0,0),1)
summary(glht(Treatment_NT, linfct = K))

# Exploratory - difference in treatment effects of news_correction and third_party correction
K <- matrix(c(0,1,-1,0,0,0,0,0,0,0,0),1)
summary(glht(Treatment_NT, linfct = K))

## Table 1 (Main Effects)---------

## Pull regression values into LaTeX (as table)
## Creates Table 1 (except for linear combinations, which are copied from lines 364-427)
screenreg(list(Treatment_BA,Treatment_NT))

## Main Effects Plot (Figure 2)--------
## Combined coefficient plot of effects of various treatments on
## belief accuracy and news outlet trust
# include only relevant outcome variables (report effects of controls textually)
predictors <- c("datafin$news_correction",
                "datafin$third_correction", "datafin$news_third_correction",
                "datafin$third_news_correction")

# order variables (instead of R default order)
name_order <- c("News outlet", "Third party", 
                "Both (news first)", "Both (third party first)",
                "Belief Accuracy Pretreatment", "News Trust Prestimulus",
                "News Trust Pretreatment", "65+", "18-24", "55-64", "Political Interest",
                "Non-Hispanic White")

# Plotting function
## Creates Figure 2
multiplot(Treatment_BA, Treatment_NT, predictors = predictors, 
          horizontal = TRUE, intercept = TRUE, by = "Coefficient", numberAngle = 0,
          interceptName = "No Correction", xlab = "Treatment effect", ylab = "Correction source", title = "",
          sort = "natural",
          newNames = setNames(c("No Correction", "News outlet", "Third party", 
                                "Both (news first)", "Both (third party first)",
                                "Belief Accuracy Pretreatment", "News Trust Prestimulus",
                                "News Trust Pretreatment", "65+", "18-24", "55-64", "Political Interest",
                                "Non-Hispanic White"), names(coefficients(Treatment_BA)))) +
  theme(text = element_text(size=10), axis.text.x = element_text(angle=0, vjust=1)) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.9, 0.9)) +
  guides(shape = "none") +
  scale_color_manual(labels = c("Belief accuracy", "News trust"), values=c("indianred2", "lightskyblue3")) +
  geom_hline(yintercept = 1.5:3.5)

## Additional results--------------
## Effect of treatments on epistemic political efficacy
Treatment_EPE <- lm(datafin$poli_efficacy ~ datafin$news_correction + datafin$third_correction
                    + datafin$news_third_correction + datafin$third_news_correction
                    + datafin$age_18 + datafin$age_25 + datafin$age_55 + datafin$age_65
                    + datafin$college_educ + datafin$nonhispanicwhite 
                    + datafin$poli_efficacy_pre)

coeftest(Treatment_EPE, vcov = vcovHC, type = "HC1")

## Regression table
## Table B2
screenreg(Treatment_EPE)

## Moderating Effects
## Effect of Republican identification on outcome measures
## Belief Accuracy x Republicans
GOP_BA <- lm(datafin$belief_accuracy ~ datafin$news_correction * datafin$Republican + datafin$third_correction * datafin$Republican
             + datafin$news_third_correction * datafin$Republican + datafin$third_news_correction * datafin$Republican
             + datafin$age_65  + datafin$belief_accuracy_pre) 
coeftest(GOP_BA, vcov = vcovHC, type = "HC1")

## News outlet trust x Republicans
GOP_NT <- lm(datafin$news_trust ~ datafin$news_correction * datafin$Republican + datafin$third_correction * datafin$Republican
             + datafin$news_third_correction * datafin$Republican + datafin$third_news_correction * datafin$Republican
             + datafin$age_18 + datafin$age_55
             + datafin$poli_interest + datafin$nonhispanicwhite 
             + datafin$news_trust_pre + datafin$news_trust_prestimulus) 
coeftest(GOP_NT, vcov = vcovHC, type = "HC1")

## Table B3
screenreg(list(GOP_BA, GOP_NT))

## Twitter usage
## Belief accuracy x Twitter usage
TW_BA <- lm(datafin$belief_accuracy ~ datafin$news_correction * datafin$twitter_user + datafin$third_correction * datafin$twitter_user
            + datafin$news_third_correction * datafin$twitter_user + datafin$third_news_correction * datafin$twitter_user
            + datafin$age_65  + datafin$belief_accuracy_pre) 
coeftest(TW_BA, vcov = vcovHC, type = "HC1")

## News outlet trust x Twitter usage
TW_NT <- lm(datafin$news_trust ~ datafin$news_correction * datafin$twitter_user + datafin$third_correction * datafin$twitter_user
            + datafin$news_third_correction * datafin$twitter_user + datafin$third_news_correction * datafin$twitter_user
            + datafin$age_18 + datafin$age_55
            + datafin$poli_interest + datafin$nonhispanicwhite 
            + datafin$news_trust_pre + datafin$news_trust_prestimulus) 
coeftest(TW_NT, vcov = vcovHC, type = "HC1")

## Export Twitter moderation regression tables
## Table B4
screenreg(list(TW_BA, TW_NT))

## General media trust moderation
# Belief accuracy x General media trust
GMT_BA <- lm(datafin$belief_accuracy ~ datafin$news_correction * datafin$media_trust + datafin$third_correction * datafin$media_trust
             + datafin$news_third_correction * datafin$media_trust + datafin$third_news_correction * datafin$media_trust
             + datafin$age_65  + datafin$belief_accuracy_pre) 
coeftest(GMT_BA, vcov = vcovHC, type = "HC1")

# News outlet trust x General media trust
GMT_NT <- lm(datafin$news_trust ~ datafin$news_correction * datafin$media_trust + datafin$third_correction * datafin$media_trust
             + datafin$news_third_correction * datafin$media_trust + datafin$third_news_correction * datafin$media_trust
             + datafin$age_18 + datafin$age_55
             + datafin$poli_interest + datafin$nonhispanicwhite 
             + datafin$news_trust_pre + datafin$news_trust_prestimulus)
coeftest(GMT_NT, vcov = vcovHC, type = "HC1")

## Regression table
## Table B5
screenreg(list(GMT_BA, GMT_NT))

## End of document------